{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "8\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy import integrate\n",
    "import scipy as sp\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "import os\n",
    "import multiprocessing\n",
    "from joblib import Parallel, delayed\n",
    "num_cores = multiprocessing.cpu_count()\n",
    "print(num_cores)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define parameter\n",
    "- 0 is defined at gate\n",
    "\n",
    "- positive is defined as moving up from gate to sample\n",
    "\n",
    "- $\\epsilon = \\sqrt{\\epsilon_{\\bot} \\epsilon_{\\parallel}}$\n",
    "\n",
    "- $\\beta = \\sqrt{\\frac{\\epsilon_{\\parallel}}{\\epsilon_{\\bot}}}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Coulomb energy (epsion=1): 137.42, magnetic length: 10.48\n",
      "beta = 1.48, epsilon = 4.45\n"
     ]
    }
   ],
   "source": [
    "#physical constant\n",
    "B_field = 6.\n",
    "EC = 56.1*np.sqrt(B_field)           #in meV (using epsilon=1, the coulomb energy in vacuum)\n",
    "lB = 25.66/np.sqrt(B_field)          #Magnetic length in nm\n",
    "\n",
    "#sample-impurity distance\n",
    "d_si = 0 / lB                  #(above: >0, below: <0)\n",
    "\n",
    "#other distance\n",
    "d_gs = 70. / lB                 #gate-sample distance\n",
    "d_sd = 2.5 / lB                 #sample-detector distance\n",
    "d_gi = d_gs + d_si              #gate to impurity distance\n",
    "d_id = d_sd - d_si              #impurity to detector distance\n",
    "\n",
    "# Dielectric constants\n",
    "beta = np.sqrt(6.6/3.)  #6.6 or 5.33\n",
    "eps  = np.sqrt(6.6*3)\n",
    "\n",
    "#smearing factor (convolve gaussian with impurity potential): sigma = lb*factor\n",
    "\n",
    "factor = 2**0.5\n",
    "\n",
    "print('Coulomb energy (epsion=1): {0:.2f}, magnetic length: {1:.2f}'.format(EC,lB))\n",
    "print('beta = {0:.2f}, epsilon = {1:.2f}'.format(beta,eps))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>filling</th>\n",
       "      <th>label</th>\n",
       "      <th>0</th>\n",
       "      <th>2.4609375</th>\n",
       "      <th>4.921875</th>\n",
       "      <th>7.3828125</th>\n",
       "      <th>9.84375</th>\n",
       "      <th>12.3046875</th>\n",
       "      <th>14.765625</th>\n",
       "      <th>17.2265625</th>\n",
       "      <th>...</th>\n",
       "      <th>100.8984375</th>\n",
       "      <th>103.359375</th>\n",
       "      <th>105.8203125</th>\n",
       "      <th>108.28125</th>\n",
       "      <th>110.7421875</th>\n",
       "      <th>113.203125</th>\n",
       "      <th>115.6640625</th>\n",
       "      <th>118.125</th>\n",
       "      <th>120.5859375</th>\n",
       "      <th>123.046875</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>-0.5</td>\n",
       "      <td>-1/2</td>\n",
       "      <td>-0.110195</td>\n",
       "      <td>-0.110657</td>\n",
       "      <td>-0.110195</td>\n",
       "      <td>-0.110506</td>\n",
       "      <td>-0.110364</td>\n",
       "      <td>-0.11008</td>\n",
       "      <td>-0.109899</td>\n",
       "      <td>-0.109502</td>\n",
       "      <td>...</td>\n",
       "      <td>-0.106549</td>\n",
       "      <td>-0.106545</td>\n",
       "      <td>-0.106516</td>\n",
       "      <td>-0.106517</td>\n",
       "      <td>-0.106513</td>\n",
       "      <td>-0.106518</td>\n",
       "      <td>-0.106505</td>\n",
       "      <td>-0.106524</td>\n",
       "      <td>-0.106525</td>\n",
       "      <td>-0.106539</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>1 rows × 53 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   filling label         0  2.4609375  4.921875  7.3828125   9.84375  \\\n",
       "0     -0.5  -1/2 -0.110195  -0.110657 -0.110195  -0.110506 -0.110364   \n",
       "\n",
       "   12.3046875  14.765625  17.2265625  ...  100.8984375  103.359375  \\\n",
       "0    -0.11008  -0.109899   -0.109502  ...    -0.106549   -0.106545   \n",
       "\n",
       "   105.8203125  108.28125  110.7421875  113.203125  115.6640625   118.125  \\\n",
       "0    -0.106516  -0.106517    -0.106513   -0.106518    -0.106505 -0.106524   \n",
       "\n",
       "   120.5859375  123.046875  \n",
       "0    -0.106525   -0.106539  \n",
       "\n",
       "[1 rows x 53 columns]"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# This data is generated by Fig2_c_1.m\n",
    "# The matlab code will save the raw data into workspace\n",
    "# Save the workspace data into an CSV file and you will get 'one_half_line.csv'\n",
    "df = pd.read_csv(r'one_half_line.csv')\n",
    "\n",
    "# remove data too far from the core\n",
    "r_range = 125\n",
    "to_pop = [c for c in df.columns if (c.replace('.','0').isnumeric() and float(c)>r_range) ]\n",
    "for c in to_pop:\n",
    "    df.pop(c)\n",
    "\n",
    "nu = np.asarray(df.iloc[:,0])                       #np array of filling factor\n",
    "label = list(df['label'])                           #np array of string of filling factor\n",
    "\n",
    "r = df.columns[2:].to_numpy().astype(float) / lB    #np array in magnetic length\n",
    "r_range = r_range / lB\n",
    "\n",
    "n_r , n_nu = (r.shape[0],nu.shape[0])               #number of r point, filling factor\n",
    "n_q = n_r*2\n",
    "\n",
    "q = np.linspace(0,10,n_q)\n",
    "q_range = q.max()\n",
    "\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Work flow:\n",
    "\n",
    "1. $\\phi_{total}(q)$ , $\\phi_{total} (r)$ \n",
    "    - Get $\\phi_{total} (r)$ from data\n",
    "\n",
    "    - $\\phi_{total}(q) = 2\\pi \\int_{0}^{\\infty}{\\phi_{total}(r)J_0(qr)rdr}$\n",
    "    \n",
    "\n",
    "2. $\\phi_{impurity}(q)$ , $\\phi_{impurity}(r)$\n",
    "\n",
    "    - $n_{impurity}(q) = 1$ $( n_{impurity}(r) = \\frac{1}{2\\pi r}\\delta(r) )$\n",
    "\n",
    "    - $\\phi_{impurity}(q) = n_{impurity}(q) \\cdot G(q) $\n",
    "\n",
    "    - $\\phi_{impurity}(r) = \\frac{1}{2\\pi} \\int_{0}^{\\infty}{n_{impurity}(q) \\cdot G(q) \\cdot J_0(qr) \\cdot N(0,factor\\times l_B) q dq}$\n",
    "\n",
    "3. $\\phi_{2DEG}(q)$ , $\\phi_{2DEG}(r)$\n",
    "\n",
    "    - $\\phi_{2DEG}(q) = \\phi_{total}(q) - \\phi_{impurity}(q)$\n",
    "    \n",
    "4. $n_{2DEG}(r)$ , $n_{2DEG}(q)$\n",
    "\n",
    "    - $n_{2DEG}(q) = \\phi_{2DEG}(q) / G(q)$\n",
    "    \n",
    "    - $n_{2DEG}(r) = \\frac{1}{2\\pi} \\int_{0}^{\\infty}{\\frac{\\phi_{2DEG}(q)}{G(q)}J_0(qr)qdq}$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Green's function\n",
    "- $G(q) = \\frac{2\\pi}{q} \\frac{ sinh(\\beta qd_g)csch(\\beta q(d_g+d_s)) }{1 + \\epsilon coth(\\beta q (d_g + d_s))} $\n",
    "- $d_g$: distance between sample/impurity and gate\n",
    "- $d_s$: distance between sample/impurity and detector\n",
    "\n",
    "    *$G(q=0)=\\frac{d_g\\beta}{\\epsilon}$ $\\newline$\n",
    "    *2D free space: $G(q) = \\frac{2\\pi}{q}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "#G(q)\n",
    "green = lambda dg,ds,iq: dg*beta/eps*2*np.pi if iq == 0. \\\n",
    "             else 2*np.pi/iq * ( np.exp(-1*iq*beta*(ds)) - np.exp( -1*iq*beta*(2*dg+ds) ) ) / (1 + eps - (1-eps) * np.exp( -2*iq*beta*(dg+ds) ))\n",
    "\n",
    "#q * G(q)\n",
    "q_green = lambda dg,ds,iq: 0 if iq == 0. \\\n",
    "             else 2*np.pi*( np.exp(-1*iq*beta*(ds)) - np.exp( -1*iq*beta*(2*dg+ds) ) ) / (1 + eps - (1-eps) * np.exp( -2*iq*beta*(dg+ds) ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "def J0_integral(rV, q, tol=1e-15, d = 1, dr = np.inf, full_output=False):\n",
    "\t\t\"\"\"MPZ: Integrates  Vq = 2 \\\\pi \\\\int_0^\\\\inf  J_0(q r) rV(r) dr.\n",
    "\t\t\t\trV is a callable function, with domain [0, +inf) and range in float.\n",
    "\t\t\t\tq is a float\n",
    "\t\t\t\ttol = Desired absolute error\n",
    "\t\t\t\td = Length scale at which any r-->0 singularities may occur\n",
    "\t\t\t\tdr = Optional; length scale beyond which rV(r) is negligible. Can be important to know when 1 >> dr q\n",
    "\t\t\t\tfull_output is a Boolean: then return number of periods integrated and #func evaluations\n",
    "\n",
    "\t\t\tReturns\n",
    "\t\t\t\tif full_output is False:\n",
    "\t\t\t\t\tVq\n",
    "\t\t\t\tif full_output is True:\n",
    "\t\t\t\t\tVq, number of periods integrated, number of function evaluation\n",
    "\n",
    "\t\t\tExample:\n",
    "\t\t\t\t>>>\tq = 2.\n",
    "\t\t\t\t>>>\tprint mod.mQHfunc.J0_integral(lambda r: np.exp(-r), q)\n",
    "\t\t\t\t>>>\tprint 2 * np.pi / np.sqrt(q**2 + 1)\n",
    "\t\t\t\t\t2.80992589242\n",
    "\t\t\t\t\t2.80992589242\n",
    "\t\n",
    "\t\t\tStrategy:\n",
    "\t\t\t1)\tIntegrate in subintervals [r_j, r_{j+1}], defined by r_j = pi (j + 3/4) / q .\n",
    "\t\t\t\t\tThe r_j's are the approximate zeros of J_0(q r),\n",
    "\t\t\t2)\tTabulate the partial sums, and Euler sum them to accelerate convergence.\n",
    "\t\t\t\t\n",
    "\t\t\t\t\n",
    "\t\t\"\"\"\n",
    "\t\ttwo_pi = 2. * np.pi\n",
    "\t\t\n",
    "\t\tdef euler_sum(S):\n",
    "\t\t\t#Euler summation given partial sums S\n",
    "\t\t\tN = len(S) - 1\n",
    "\t\t\tbinom = [ sp.special.comb(N, n, exact = True) for n in range(N+1) ]\t\t# binomial coef\n",
    "\t\t\treturn np.inner(S, binom) / 2**N\n",
    "\t\t\t#return np.inner(S, sp.stats.binom.pmf(range(N+1),N, 0.5))\n",
    "\t\tdef I(r):\n",
    "\t\t\t#The integrand\n",
    "\t\t\treturn sp.special.j0(q*r)*rV(r)\n",
    "\t\n",
    "\t\tn_Eval = 0\t\t# number of rV() evaluations\n",
    "\t\tn_Quad = 0\n",
    "\t\tstatus = 'converged'\n",
    "\t\tif q==0.:\n",
    "\t\t\t#Break out near singular behaviour near 0\n",
    "\t\t\ta0, abserr, infodict = sp.integrate.quad(rV, 0., d, epsabs=tol/50., full_output = True)[:3]\n",
    "\t\t\tn_Eval+=infodict['neval']\n",
    "\n",
    "\t\t\ta1, abserr, infodict = sp.integrate.quad(rV, d, np.inf, epsabs=tol/50., full_output = True)[:3]\n",
    "\t\t\tn_Eval+=infodict['neval']\n",
    "\t\t\tn_Quad+=2\n",
    "\t\t\tif full_output:\n",
    "\t\t\t\treturn (a0+a1) * two_pi, {'n_quad':n_Quad, 'n_eval':n_eval, 'converged':True}\n",
    "\t\t\telse:\n",
    "\t\t\t\treturn (a0+a1) * two_pi\n",
    "\t\t\t\n",
    "\t\t#Seperate out first interval, [0, r_0], and garauntee subdivision for accuracy\n",
    "\t\tbreak_pt = min(0.375*np.pi/q, d)\n",
    "\t\tif np.pi*0.75/q < dr:\n",
    "\t\t\tbreak_pt = min(0.375*np.pi/q, d)\n",
    "\t\t\ta0, abserr, infodict =  sp.integrate.quad(I, 0., np.pi*0.75/q, epsabs=tol/50., full_output = True, points = [break_pt])[:3]\n",
    "\t\t\tn_Eval+=infodict['neval']\n",
    "\t\t\tn_Quad+=1\n",
    "\t\telse:\n",
    "\t\t\tbreak_pt = min(dr/2., d)\n",
    "\t\t\ta0, abserr, infodict =  sp.integrate.quad(I, 0., dr, epsabs=tol/50., full_output = True, points = [break_pt])[:3]\n",
    "\t\t\tn_Eval+=infodict['neval']\n",
    "\t\t\ta1, abserr, infodict =  sp.integrate.quad(I, dr, np.pi*0.75/q, epsabs=tol/50., full_output = True)[:3]\n",
    "\t\t\tn_Eval+=infodict['neval']\n",
    "\t\t\tn_Quad+=2\n",
    "\t\t\ta0+=a1\n",
    "\t\ta0 *= two_pi\n",
    "\t\t\n",
    "\t\t#Initialize S, ES\n",
    "\t\tr, abserr, infodict = sp.integrate.quad(I, np.pi*0.75/q, np.pi*1.75/q, full_output = True)[:3]\n",
    "\t\tr *= two_pi\n",
    "\t\tn_Eval+=infodict['neval']\n",
    "\t\tn_Quad+=1\n",
    "\t\tS = [r] #partial sums\n",
    "\t\tES = [r] #euler sums\n",
    "\t\t\n",
    "\t\tj = 1.75\n",
    "\t\twhile j < 50:\n",
    "\t\t\tr, abserr, infodict = sp.integrate.quad(I, np.pi*j/q, np.pi*(j+1.)/q, full_output = True)[:3]\n",
    "\t\t\tr *= two_pi\n",
    "\t\t\tn_Eval+=infodict['neval']\n",
    "\t\t\tn_Quad+=1\n",
    "\t\t\tS.append( S[-1] + r )\n",
    "\t\t\t\n",
    "\t\t\tif abs(r) < max(abserr, tol): #Forget Euler summation - absolute convergence!\n",
    "\t\t\t\tif full_output:\n",
    "\t\t\t\t\treturn S[-1] + a0, {'n_quad':n_Quad, 'n_eval':n_eval, 'converged':True}\n",
    "\t\t\t\telse:\n",
    "\t\t\t\t\treturn S[-1] + a0\n",
    "\t\t\tES.append(euler_sum(S))\n",
    "\t\t\n",
    "\t\t\tif\tabs(ES[-1] - ES[-2]) < max(abserr, tol): #Good enough - our work is done.\n",
    "\t\t\t\tif full_output:\n",
    "\t\t\t\t\treturn ES[-1] + a0, {'n_quad':n_Quad, 'n_eval':n_eval, 'converged':True}\n",
    "\t\t\t\treturn ES[-1] + a0\n",
    "\t\t\t\n",
    "\t\t\tj+=1.\n",
    "\t\n",
    "\t\tprint(\"\\tJ0_integral: Poorly behaved integrand!  q = {}, d = {}, tol = {}, n_Quad = {}, n_Eval = {}, \".format(q, d, tol, n_Quad, n_Eval))\n",
    "\t\tprint(\"\\t\\tS - S[-1]:\", np.array(S) - S[-1])\n",
    "\t\tprint(\"\\t\\tES - ES[-1]:\", np.array(ES) - ES[-1])\n",
    "\t\tif full_output:\n",
    "\t\t\treturn ES[-1] + a0, {'n_quad':n_Quad, 'n_eval':n_eval, 'converged':False, 'dS':np.array(S) - S[-1], 'dES':np.array(ES) - ES[-1]}\n",
    "\t\telse:\n",
    "\t\t\treturn ES[-1] + a0\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "#create callalbe function for J0_integral\n",
    "class r_intep():\n",
    "    def __init__(self,data):\n",
    "        self.tmp_f = data\n",
    "    def rf(self,r_query):\n",
    "        if r_query > r_range: return 0\n",
    "        else: return r_query*np.interp(r_query , r.tolist() , self.tmp_f.tolist())\n",
    "\n",
    "class q_intep():\n",
    "    def __init__(self,data):\n",
    "        self.tmp_f = data\n",
    "    def qf(self,q_query):\n",
    "        if q_query > q_range: return 0\n",
    "        else: return q_query*np.interp(q_query , q.tolist() , self.tmp_f.tolist())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. $\\phi_{total}(q)$ , $\\phi_{total} (r)$ \n",
    "    \n",
    "- Get $\\phi_{total} (r)$ from data\n",
    "\n",
    "- $\\phi_{total}(q) = 2\\pi \\int_{0}^{\\infty}{\\phi_{total}(r)J_0(qr)rdr}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "number of data point (r):  51\n",
      "number of data in each LL branch:  1\n"
     ]
    }
   ],
   "source": [
    "phi_r_total = []\n",
    "phi_r_total_func = []\n",
    "mu = []\n",
    "\n",
    "#\\phi_{impurity}(r)\n",
    "for i in range(n_nu):\n",
    "    tmp = df.iloc[i,2:].to_numpy()*1000 / EC\n",
    "    mu.append(np.mean(tmp[-10:]))\n",
    "    tmp = tmp-mu[i]\n",
    "    phi_r_total.append(tmp)         #.reshape(-1, 1)\n",
    "\n",
    "#\\phi_{impurity}(q)\n",
    "phi_r_total_func = [r_intep(phi_r_total[id]) for id in range(n_nu)]\n",
    "\n",
    "run_phi_q_total = lambda ind: np.array([J0_integral(phi_r_total_func[ind].rf, iq,\\\n",
    "                                                    tol=1e-21, d = 1, dr = r_range,\\\n",
    "                                                    full_output=False) for iq in q])\n",
    "\n",
    "phi_q_total = Parallel(n_jobs=num_cores)(delayed(run_phi_q_total)(i) for i in range(n_nu))\n",
    "\n",
    "print('number of data point (r): ',n_r)\n",
    "print('number of data in each LL branch: ',n_nu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. $\\phi_{impurity}(q)$ , $\\phi_{impurity}(r)$\n",
    "- $n_{impurity}(r) = \\frac{1}{2\\pi r}\\delta(r)$\n",
    "\n",
    "- $n_{impurity}(q) = 1$\n",
    "\n",
    "- $\\phi_{impurity}(q) = n_{impurity}(q) \\cdot G(q) $\n",
    "\n",
    "\n",
    "- $\\phi_{impurity}(r) = \\frac{1}{2\\pi} \\int_{0}^{\\infty}{n_{impurity}(q) \\cdot G(q) \\cdot J_0(qr)\\cdot N(0,factor\\times l_B)qdq}$\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\\n_{impurity}(q)\n",
    "n_q_imp = -1\n",
    "\n",
    "#G(q)_{impurity}\n",
    "green_imp = lambda iq: green(d_gi,d_id,iq)\n",
    "\n",
    "#q* G(q)_{impurity}\n",
    "q_green_imp = lambda iq: q_green(d_gi,d_id,iq)\n",
    "\n",
    "#\\phi_{impurity}(q)\n",
    "phi_q_imp = lambda iq: n_q_imp * green_imp(iq)\n",
    "\n",
    "#q * \\phi_{impurity}(q)\n",
    "q_phi_q_imp = lambda iq: n_q_imp * q_green_imp(iq)\n",
    "\n",
    "#q * convolve( \\phi_{impurity}(q) )\n",
    "q_phi_q_imp_convolve = lambda iq: q_phi_q_imp(iq) * np.exp(-0.5*(iq*factor)**2) \n",
    "\n",
    "#convolve( \\phi_{impurity}(q) )\n",
    "phi_q_imp_convolve = lambda iq: phi_q_imp(iq) * np.exp(-0.5*(iq*factor)**2)\n",
    "\n",
    "# \\phi_impurity(r)\n",
    "phi_r_imp = np.array([J0_integral(q_phi_q_imp_convolve, ir,tol=1e-21, d = 1,\\\n",
    "                                  dr = np.inf,full_output=False) for ir in r]) / (2*np.pi)**(2)\n",
    "\n",
    "\n",
    "df_green = pd.DataFrame()\n",
    "df_green['q'] = q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. $\\phi_{2DEG}(q)$ , $\\phi_{2DEG}(r)$\n",
    "\n",
    "- $\\phi_{2DEG}(q) = \\phi_{total}(q) - \\phi_{impurity}(q)$\n",
    "\n",
    "- $\\phi_{2DEG}(r) = \\phi_{total}(r) - \\phi_{impurity}(r)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "phi_q_2DEG = [phi_q_total[id] - [phi_q_imp_convolve(iq) for iq in q] for id in range(n_nu)]\n",
    "phi_r_2DEG = [phi_r_total[id] - phi_r_imp for id in range(n_nu)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. $n_{2DEG}(r)$ , $n_{2DEG}(q)$\n",
    "\n",
    "- $n_{2DEG}(q) = \\phi_{2DEG}(q) / G(q)$\n",
    "\n",
    "- filter out high q signal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "green_sample = np.array([green(d_gs,d_sd,iq) for iq in q])\n",
    "\n",
    "#filter\n",
    "z,p,k = sp.signal.butter(4, 3/max(q), 'low', output='zpk' )\n",
    "w, h = sp.signal.freqz_zpk(z, p, k , worN=n_q)\n",
    "h = abs(h)\n",
    "\n",
    "n_q_2DEG = [phi_q_2DEG[id] / green_sample * h for id in range(n_nu)]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- $n_{2DEG}(r) = \\frac{1}{2\\pi} \\int_{0}^{\\infty}{\\frac{\\phi_{2DEG}(q)}{G(q)}J_0(qr)qdq}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_q_2DEG_func = [q_intep(n_q_2DEG[id]) for id in range(n_nu)]\n",
    "\n",
    "run_n_r_2DEG = lambda ind: np.array([J0_integral(n_q_2DEG_func[ind].qf, ir,\\\n",
    "                                                    tol=1e-21, d = 1, dr = r_range,\\\n",
    "                                                    full_output=False) for ir in r]) / (2*np.pi)**(2)\n",
    "\n",
    "n_r_2DEG = Parallel(n_jobs=num_cores)(delayed(run_n_r_2DEG)(i) for i in range(n_nu))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_save = pd.DataFrame()\n",
    "df_save['position'] = r\n",
    "df_save['phi_r_total'] = phi_r_total[0]\n",
    "df_save['phi_r_imp'] = phi_r_imp\n",
    "df_save['phi_r_2DEG'] = phi_r_2DEG[0]\n",
    "df_save['n_r_2DEG'] = n_r_2DEG[0]\n",
    "df_save.to_csv(r'result_one_half.csv')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Princeton",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
